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Equilibrium  and  non-equilibrium  molecular-dynamics  simulations  are  employed  in  this  study  to  investigate 
various  aspects  of  shock  waves  in  fused  silica  (a  pure  Si02  amorphous  material  used  in  transparent-armor 
applications).  Equilibrium  molecular-dynamics  simulations  are  used  first  to  validate  that  the  initial  (un¬ 
shocked)  fused  silica  possesses  the  appropriate  mass  density  and  microstructure  (as  characterized  by  its 
partial  Si-Si,  Si-O,  and  0-0  radial  distribution  functions).  Next,  non-equilibrium  molecular-dynamics 
simulations  are  employed,  within  a  continuously  contracting  computational-cell  scheme,  to  generate  planar 
longitudinal  (uniaxial  motion)  shocks  of  different  strengths.  By  examining  and  quantifying  the  dynamics  of 
shock-wave  motion,  the  respective  shock-Hugoniot  relations  (i.e.,  functional  relations  between  various 
material-state  variables  in  the  material  states  produced  by  the  shocks  of  different  strengths)  are  determined. 
This  methodology  suggested  that  irreversible  non-equilibrium  deformation/damage  processes  play  an 
important  role  in  the  mechanical  response  of  fused  silica  to  shock  loading  and  that  the  “equilibrium” 
procedures  for  Hugoniot  determination  based  on  the  equation  of  state  and  the  Rankine-Hugoniot  equation 
may  not  be  fully  justified.  Finally,  the  non-equilibrium  molecular-dynamics  simulations  were  used  to 
identify  the  main  microstructure  modifying/altering  processes  accompanying  the  shock-wave  motion 
through  fused  silica. 
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1.  Introduction 

Transparent  blast/ballistic-impact-resistant  vehicle  structures 
(e.g.,  windshields,  door  windows,  viewports,  etc.),  employed  in 
ground  and  air  vehicles,  utilize  a  variety  of  non-glass  transpar¬ 
ent  materials  such  as  transparent  crystalline  ceramics  (e.g., 
aluminum-oxinitride  spinel,  A10N,  sapphire  (Ref  1)),  and  new 
transparent  polymer  materials  (e.g.,  transparent  nylon  (Ref  2)). 
Despite  the  well-established  benefits  offered  by  these  novel 
transparent  materials,  ballistic  ceramic-based  glasses  are  still 
important  constituent  materials  in  a  majority  of  transparent 
impact-resistant  structures  (i.e.,  transparent  armor)  used  today. 
The  main  reasons  for  the  wide  usage  of  ceramic  glasses  are  (a) 
glass- structure  fabrication  technologies  enable  the  production 
of  curved,  large  surface-area,  transparent  structures  with 
thickness  approaching  several  inches;  (b)  material  processing 
and  product  manufacturing  are  associated  with  relatively  low 
costs;  and  (c)  ballistic  and  blast  impact  survivability  of  ceramic 
glasses  can  be  substantially  improved  through  the  use  of 
techniques,  such  as  compositional  modifications,  chemical 
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strengthening,  controlled  crystallization,  and  through  the  use  of 
various  fabrication-process  control  strategies  (Ref  3). 

In  addition  to  their  good  overall  blast/ballistic-impact 
capabilities,  some  classes  of  ceramic  glasses  have  been  found 
to  exhibit  particularly  good  performance  against  shaped  charge 
jets.  These  (molten-metal,  high-velocity)  jets  are  produced 
during  impact  of  warheads  consisting  of  an  explosive  (located 
at  the  back  of  the  warhead  and  capped  at  the  front  by  a 
backward  pointing  metal  liner)  and  a  conical  hollow  front 
section  (Ref  4).  The  extremely  high  explosion-induced  pres¬ 
sures  cause  the  forward  inversion  of  the  metal  liner  along  the 
warhead  axis,  and  the  accompanying  large  dissipative  plastic 
work  results  in  the  formation  of  a  high-velocity  forward 
traveling  jet(s)  of  molten  metal.  While  the  aforementioned 
finding  pertaining  to  the  superior  resistance  of  ceramic  glasses 
to  impact  by  shaped  charge  jets  are  believed  to  be  closely 
related  to  the  glass  chemistry  and  micro  structure  of  these 
materials,  details  of  the  underlying  chemistry/micro  structure/ 
performance  relations  are  not  well  understood. 

The  development  of  new  glass-based,  transparent,  impact- 
resistant  structures,  aimed  at  reducing  the  vulnerability  of 
protected-vehicle  occupants  and  on-board  instrumentation  to 
various  blast/ballistic  threats,  is  nowadays  mainly  based  on  the 
empiricism,  legacy  knowledge,  and  the  traditional  “fabricate- 
and-test”  strategies.  While  these  experimental  strategies  are 
very  critical  for  ensuring  the  reliability  and  effectiveness  of  the 
transparent,  impact-resistant  structures,  they  are  generally 
associated  with  high  costs,  long  lead  times  and  destructive 
“one-shot”  testing.  In  recent  years,  these  experimental  strate¬ 
gies  have  begun  to  be  increasingly  complemented  by  the 
corresponding  computation-based  modeling  and  simulation 
efforts.  However,  it  is  now  well-established  (Ref  5-7)  that  the 
effectiveness  and  reliability  of  the  computation-based  modeling 
and  simulation  approaches  are  greatly  affected  by  the  fidelity  of 


Journal  of  Materials  Engineering  and  Performance 


Volume  21(6)  June  2012—823 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

JUN  2012  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2012  to  00-00-2012 

4.  TITLE  AND  SUBTITLE 

Molecular-Level  Analysis  of  Shock-Wave  Physics  and  Derivation  of  the 
Hugoniot  Relations  for  Fused  Silica 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS (ES) 

Clemson  University , Department  of  Mechanical  Engineering, 241 
Engineering  Innovation  Building, Clemson, SC, 29634 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS  (ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

ARSTRATT 

1 8 .  NUMBER  1 9a.  NAME  OF 

OF  PAGES  RESPONSIBLE  PERSON 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  Same  aS 

unclassified  unclassified  unclassified  Report  (SAR) 

14 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


the  employed  material  models  (mathematical  relations  em¬ 
ployed  to  represent  interdependencies  among  various  material 
state  variables).  In  the  case  of  transparent  protection  systems 
under  investigation  in  this  study,  it  is  critical  that  these  models 
realistically  describe  deformation/fracture  response  of  the 
subject  material  (ballistic  (fused  silica-based)  glass,  in  the 
present  case)  under  high-rate,  large-strain,  and  high-pressure 
loading  conditions  encountered  during  blast/ballistic  impact. 
Therefore,  one  of  the  main  objectives  of  the  present  study  is  to 
further  advance  the  application  of  computational  modeling/ 
simulation-based  engineering  approaches  of  transparent  im¬ 
pact-resistant  structures  via  the  potential  improvements  in 
fidelity  of  the  associated  material  (fused  silica)  models  (which, 
in  turn,  will  be  achieved  via  identification  and  quantification  of 
processes  and  phenomena  occurring  in  fused  silica  under  such 
conditions.  Specifically,  phenomena  and  processes  associated 
with  the  formation  and  propagation  of  shock-waves  within 
fused  silica  will  be  investigated. 

A  comprehensive  literature  review  carried  out  as  part  of  our 
prior  studies  (Ref  8,  9)  revealed  that  the  mechanical  behavior  of 
glass  has  been  modeled  predominantly  using  three  distinct 
approaches:  (a)  molecular-modeling  methods  (Ref  8-14);  (b) 
continuum-material  approximations  (Ref  5-7,  15-18);  and  (c) 
meso-length  scale  models  based  on  explicit  representation  of 
cracks  (Ref  19,  20).  Since  a  detailed  overview  of  these  models 
can  be  found  in  Ref  8  and  9,  they  will  not  be  discussed  any 
further  in  the  present  article.  However,  it  is  worth  mentioning 
that  the  overviews  presented  in  Ref  8  and  9  clearly  established 
that  (a)  molecular-level  models  are  critical  for  identifying 
nanometer  length-scale  phenomena  and  the  associated  micro¬ 
structure  evolution  processes  accompanying  shock-wave  prop¬ 
agation  and  damage/failure  evolution;  (b)  continuum-level 
models  are  the  only  ones  which  enable  large-scale  computa¬ 
tional  analyses  of  the  mechanical  response  of  full-scale 
transparent-armor  protective  structures  to  various  blast  and 
ballistic  threats;  and  (c)  the  models  based  on  explicit  represen¬ 
tation  of  cracks  and  their  propagation  are  mainly  suitable  for 
the  study  of  the  local  deformation  and  fracture  response  of 
simple-geometry  test  coupons. 

Ceramic  glasses  (like  fused  silica,  investigated  in  the  present 
study)  are  amorphous  materials.  The  molecular-level  micro¬ 
structure  of  these  materials  involves  entities,  such  as  a  random 
network  of  covalently  bonded  atoms,  atomic  free  volumes, 
(network-former)  bridging  and  non-bridging  oxygen  atoms, 
cations  of  glass-modifier  species,  etc.  Despite  the  absence  of  a 
crystalline  structure,  the  microstructure  of  ceramic  glasses  is 
not  completely  random  but  rather  involves  different  extents  of 
short-  and  intermediate-range  order,  spanning  over  a  range  of 
length-scales  (from  the  quantum-mechanical  to  the  continuum- 
level).  That  is  the  reason  that  as  a  part  of  this  larger  research 
effort,  involving  the  scientists  from  the  Army  Research 
Laboratory,  Aberdeen  Proving  Ground,  MD  and  Clemson 
University,  a  multi  length-scale  modeling  and  simulation 
approaches  have  been  adopted  in  the  investigation  of  various 
ceramic  glasses. 

A  schematic  of  the  adopted  multi-length-scale  approach 
is  depicted  in  Fig.  1  in  which  the  characteristic  length-range  is 
shown  along  the  abscissa,  and  the  characteristic  time-scale  is 
shown  along  the  ordinate  for  the  following  selected  length- 
scales:  (a)  electronic  scale:  within  which  the  material  is 
investigated  using  quantum-mechanical  methods  and  tools; 
(b)  atomic/molecular  length-scale — within  which  molecular- 
statics  and  -dynamics  computational  techniques  are  employed; 
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Fig.  1  A  schematic  of  the  multi-scale  modeling  and  simulation 
approach  utilized  in  our  ongoing  study  on  ceramic  glasses 


(c)  micro  length- scale — within  which  various  aspects  of  the 
short-range  order  are  investigated  using  Monte  Carlo-based 
stochastic  methods;  (d)  meso  length-scale — within  which  the 
role  of  intermediate-range  order  and  long-range  disorder  is 
analyzed  using  various  network-dynamics  computational 
schemes;  and  (e)  continuum  length-scale — within  which  the 
contribution  of  all  the  finer  length-scales  is  incorporated  into 
the  appropriate  homogenized/smeared-out  material  model,  and 
subsequently  used  within  finite  element-based  transient  non¬ 
linear  dynamics  computational  methods  and  tools.  While  all  the 
aforementioned  length-scales  are  being  investigated  in  our 
ongoing  research  and  several  classes/grades  of  ceramic  glasses 
are  being  analyzed,  the  present  article  deals  only  with  the 
molecular-level  length-scale  and  with  fused  silica. 

As  mentioned  above,  within  the  present  study,  molecular- 
level  computational  methods  are  employed  to  investigate 
shock-wave-related  phenomena  in  fused  silica.  A  shock  wave 
(or  shock  for  short)  is  a  wave  which  propagates  through  a 
medium  at  a  speed  higher  than  the  sound  speed,  and  its  passage 
causes  an  abrupt  (discontinuity-like)  change  in  the  material 
state  variables  (e.g.,  pressure/stress,  internal  energy  density, 
mass-density,  temperature,  and  particle  velocity).  Shocks  of 
higher  strength  give  rise  to  larger  extents  of  state-variable 
changes  and  propagate  at  higher  speeds.  Furthermore,  unlike 
the  acoustic  waves  which  produce  isentropic  (constant  entropy) 
changes  in  the  material  state,  shocks  bring  about  irreversible 
(entropy-increasing/energy-dissipating)  material- state  changes. 
In  fluids,  energy  dissipation  is  mainly  caused  by  high-rate 
viscous  damping  phenomena,  while  in  solids,  inelastic  defor¬ 
mation  processes  play  a  key  role  in  the  irreversibility  of 
material-state  changes  and  in  energy  dissipation  associated  with 
shock  loading. 

In  our  recent  study  (Ref  9,  21),  molecular-level  computa¬ 
tional  methods  and  tools  were  used  to  investigate  shock-related 
phenomena  in  soda-lime  glass  and  polyurea.  The  key  findings 
obtained  in  these  studies  can  be  summarized  as  follows:  (a) 
propagation  of  a  disturbance  caused  by  (impulsive)  external 
loading  typically  results  in  the  formation  of  steady  (time- 
invariant  structured)  shocks,  and  this  process  is  facilitated  by 
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the  operation  of  intrinsic  rate-dependent  energy-dissipative 
mechanisms.  However,  unlike  fluids  in  which  shock- induced 
energy  dissipation  is  mainly  related  to  viscous  damping,  in 
solids,  the  dominant  energy-dissipation  mechanisms  are  asso¬ 
ciated  with  various  inelastic-deformation/damage  processes 
resulting  from  concerted  lateral  slippage  of  the  particles;  (b)  in 
the  case  of  soda-lime  glass,  the  inelastic-deformation  processes 
involve  chemical-bond  breaking,  ring-size  modifications,  and 
atomic  coordination/network  structure  changes.  In  the  case  of 
(segmented  and  micro-segregated)  polyurea,  on  the  other  hand, 
these  processes  were  mainly  related  to  the  degradation  and 
breakage  of  the  hydrogen  bonds  within  the  so-called  hard- 
segments;  and  (c)  in  order  to  eliminate  free-surface  effects, 
molecular-level  modeling  of  shocks  had  to  be  carried  out  using 
computational  systems  with  periodic  boundary  conditions  (at 
least  in  the  directions  transverse  to  the  shock-wave  propagation 
direction).  To  attain  steady-wave  conditions  of  the  shock, 
computational  domains  sufficiently  long  in  the  direction  of 
shock-wave  propagation  had  to  be  employed.  Furthermore, 
lateral  dimensions  of  the  computational  domain  had  to  be  also 
sufficiently  large  to  avoid  spurious  effects  associated  with  the 
use  of  the  periodic  boundary  conditions.  Consequently,  com¬ 
putational  domains  involving  several  tens  of  thousands  of 
atoms  had  to  be  employed.  The  corresponding  computational 
times  (controlled  by  the  shock-wave  travel  time)  were  of  the 
order  of  5-20  ps.  This  limited  the  shock  thicknesses  to  around 
10  nm  and  the  corresponding  rise  times  to  ca.  0.5-1  ps.  Thus, 
weak  shocks  with  thicknesses  of  hundreds  of  nanometers  could 
not  be  analyzed  using  molecular-level  methods  (at  least  in  their 
steady-wave  regime). 

While  recognizing  the  aforementioned  aspects  and  potential 
limitations  of  molecular-level  modeling  of  shock,  the  present 
study  addresses  the  problem  of  shock  generation/propagation  in 
fused  silica  and  has  the  following  two  main  objectives: 

(a)  Determination  of  the  shock  Hugoniot  relations  (centered 
on  the  initial  stress-free  quiescent  state)  in  fused  silica. 
Hugoniot  is  the  (stress  versus  specific  volume  versus 
energy  density  versus  temperature  versus  entropy  density 
versus  particle  velocity  versus  shock  speed)  locus  of  the 
material  states  generated  by  the  passage  of  shocks  of 
various  strengths.  It  is  most  often  used  (i)  in  the  deriva¬ 
tion  of  the  continuum-level  high-loading  rate  material 
models  (particularly  in  the  derivation  of  the  equation  of 
state);  and  (ii)  in  the  analysis  of  planar-shock  dynamics 
and  various  shock  reflection/interaction  phenomena  and 

(b)  Characterization  and  quantification  of  the  material  state 
and  the  microstructure-altering  processes  resulting  from 
the  passage  of  shocks  of  various  strengths.  Among  these 
processes  are  those  which  involve  reversible  and  irre¬ 
versible  densifications,  chemical-bond  breaking,  coordi¬ 
nation-number  modifications,  ring-size  alterations,  and 
other  manifestations  of  the  changes  to  the  short-  and 
intermediate-range  orders.  Such  changes  have  been  ob¬ 
served  even  under  quasi-static  loading  conditions  such 
as  those  encountered  during  simple  Hertzian  indentation 
testing  (Ref  3). 

The  organization  of  the  article  is  as  follows:  A  brief 
overview  of  the  random-network  model  for  micro  structure  of 
amorphous  materials  and  its  application  to  fused  silica  are 
provided  in  section  2.  The  molecular-level  modeling  and 
simulation  procedure  employed  is  presented  in  great  detail 


in  section  3.  Key  results  pertaining  to  validation  of  the 
room-temperature  density  and  microstructure  of  fused-silica, 
dynamics  of  planar,  longitudinal  shock-wave  motion,  and 
shock-Hugoniot  relations  are  presented  and  discussed  in 
section  4.  A  summary  of  the  main  findings  and  conclusions 
is  provided  in  section  5. 


2.  Molecular-Level  Microstructure  of  Fused  Silica 

As  mentioned  earlier,  one  of  the  basic  properties  of  ceramic 
glasses  (or  of  glasses,  in  general)  is  their  lack  of  long-range 
atomic  order  which  places  them  into  a  class  of  amorphous 
materials.  For  instance,  the  atomic  arrangement  in  pure  silicate 
glass  (i.e.,  fused  silica)  is  highly  random  relative  to  the 
chemically  equivalent  (crystalline)  quartz.  To  describe  the 
structure  of  ceramic  glasses  as  determined  using  various 
experimental  techniques  (e.g.,  neutron-diffraction,  nuclear 
magnetic  resonance,  small  angle  X-ray  scattering  (SAXS), 
etc.),  the  so-called  random  network  model  (Ref  22)  is  typically 
employed.  Such  a  model  represents  amorphous  materials  as  a 
three-dimensional-linked  network  of  polyhedra.  The  character 
(number  of  facets)  of  the  polyhedra  is  controlled  by  the  species- 
specific  coordination  of  the  central  (glass-forming)  atom 
(cation).  In  the  case  of  silicate-based  glasses,  like  fused  silica 
and  soda-lime  glass,  the  polyhedron-center  atoms  are  all 
silicon,  and  each  silicon  atom  is  surrounded  by  four  oxygen 
atoms  (while  each  oxygen  atom  is  connected  to  or  bridges  two 
silicon  atoms)  forming  a  Si044-  tetrahedron.  In  glasses  of 
different  formulations,  other  types  of  polyhedra  may  exist. 
Since  silicon  has  a  tendency  to  form  a  continuous  network  with 
(bridging)  oxygen  atoms,  Si02  is  commonly  referred  to  as  a 
“network  former.”  Other  potential  network  formers  in  glass  are 
boron  and  germanium  oxides. 

Numerous  oxides  and  other  additives  are  utilized  to  modify 
the  basic  silica  tetrahedra  network  of  silicate-based  glasses  to 
tailor  their  properties  to  specific  applications.  When  alkali  (or 
alkaline  earth)  oxides  are  added  to  a  pure  silicate-based  glass,  to 
accommodate  the  excess  of  oxygen  anions  which  are  present 
due  to  the  oxide  dissociation,  the  continuity  of  the  silica 
tetrahedra  network  becomes  disrupted.  The  resulting  glass 
structure  contains  additional  non-bridging  (connected  to  only 
one  silicon  atom)  oxygen  atoms.  Charge  transfer  from  the  alkali 
earth  metal  atoms  converts  the  non-bridging  oxygen  atoms  into 
singly  charged  anions.  The  metallic  cations  formed  in  this 
process  tend  to  hover  around  the  non-bridging  oxygen  ions  for 
local  charge  neutrality.  Since  alkali  (or  alkaline  earth)-based 
oxides  cause  a  disruption  in  the  continuous  glass  network,  they 
are  typically  referred  to  as  “network  modifiers.”  As  mentioned 
earlier,  fused-silica  ceramic  glass  investigated  in  the  present 
study  is  chemically  pure  Si02  (i.e.,  free  of  glass  modifiers)  and 
contains  a  continuous  network  of  Si-0  bonds  (i.e.,  it  is  free  of 
non-bridging  oxygen  atoms). 

Within  the  random  network  model,  the  micro  structure  of 
glass  is  described  using  several  network  state  parameters. 
Among  these,  the  most  frequently  used  are  (a)  the  so-called  R 
parameter,  defined  as  the  average  number  of  oxygen  ions  per 
network-forming  ion  which  describes  the  overall  connectivity 
of  a  given  network.  In  the  case  of  fused  silica,  in  which  there 
are  two  (bridging)  oxygen  anions  for  every  network-forming 
silicon  cation,  the  R  value  is  2.0.  In  general,  a  larger  value  of 
the  R  parameter  (brought  about  by  the  addition  of  network 
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modifiers)  implies  a  more  open  (less  connected)  weaker 
amorphous  glass  network.  The  effect  of  network  formers  on 
the  R  parameter  is  more  complicated  and  depends  on  the 
network  former  coordination  number,  the  strength  of  its  bond 
with  oxygen  as  well  as  its  concentration;  (b)  the  so-called  X 
parameter  which  defines  the  average  number  of  non-bridging 
(connected  to  only  a  single  network-forming  cation)  oxygen 
atoms  per  network  polyhedron  and  takes  on  a  zero  value  in  the 
case  of  fused  silica  and;  (c)  the  so-called  Y  parameter  which 
defines  the  average  number  of  bridging  (connected  to  two 
network-forming  cations)  oxygen  atoms  per  network  polyhe¬ 
dron  and  takes  on  a  value  of  4.0  in  the  case  of  fused  silica. 

In  addition  to  the  three  parameters  mentioned  above,  the 
“seemingly”  random  micro  structure  of  ceramic  glasses  is  also 
described  using  bond-length,  bond-angle,  pair-wise  correlation, 
Voronoi-cell  volume,  etc.  distribution  functions.  A  more 
detailed  description  of  these  micro  structural  parameters  will 
be  provided  in  section  4,  as  part  of  a  discussion  of  the  results 
pertaining  to  the  structure  of  the  as-shocked  fused  silica. 

3.  Molecular-Level  Analysis  of  Fused  silica 

As  mentioned  earlier,  molecular-level  computational  meth¬ 
ods  have  been  employed  in  the  present  study  to  investigate 
various  shock-wave-related  phenomena  in  fused  silica.  Within 
these  methods,  all  atoms  and  bonds  are  explicitly  accounted  for 
and  molecular  mechanics  and  dynamics  algorithms  are 
employed  to  quantify  the  state  and  behavior  of  the  material 
under  investigation.  Molecular-level  simulation  problems  typ¬ 
ically  require  the  specification  of  the  following:  (a)  a  molecular- 
level  computational  model  consisting  of  atoms,  ions,  functional 
groups,  and/or  molecules;  (b)  a  set  of  force  field  functions;  and 
(c)  computational  method(s)  to  be  used  in  the  simulation.  More 
details  of  these  three  aspects  of  the  molecular-level  modeling 
and  simulation  of  fused  silica  are  provided  later. 

3.1  Computational  Model 

At  the  molecular  level,  fused  silica  is  modeled  as  a  discrete- 
particle-based  material  consisting  of  silicon  (Si)  and  oxygen 
(O)  atoms  mutually  bonded  via  a  single  covalent  bond  and 
forming  a  connected,  non-structured/amorphous  network  of 
silica  (Si044-)  tetrahedra. 

While  fused  silica  is  an  amorphous  material  and  does  not 
possess  any  long-range  regularity  in  its  atomic/molecular 
structure,  modeling  of  bulk  behavior  of  fused  silica  is  typically 
done  at  the  molecular  level  by  assuming  the  existence  of  a 
larger  (amorphous)  unit  cell.  Repetition  of  this  cell  in  the  three 
orthogonal  directions  (the  process  also  known  as  application  of 
the  “periodic  boundary  conditions”)  results  in  the  formation  of 
an  infinitely  large  bulk-type  material.  This  procedure  has  been 
adopted  in  the  present  study. 

The  parallelepiped- shaped  computational  cell  used  in  the 
present  (shock-wave)  analysis  contained  4608  particles  (1536  Si 
atoms  and  3072  O  atoms).  The  edge-lengths  of  the  computational 
cell  were  initially  set  to  a  =  10.36  nm  and  b  =  c  =  2.59  nm, 
yielding  a  fused  silica’s  initial  nominal  density  of  2.2  g/cm3.  The 
three  edges  (a,  b ,  and  c)  of  the  cell  were  aligned,  respectively, 
with  the  three  coordinate  axes  ( x ,  y,  and  z). 

To  create  the  ambient  temperature/pressure  equilibrium 
atomic  configuration  within  the  computational  cell,  the  follow¬ 
ing  procedure  was  implemented  within  the  Visualizer  (Ref  23) 


program  from  Accelrys:  (a)  a  starting  computational  cell  was 
first  constructed  by  stacking  the  appropriate  number  of  fused 
silica  unit  cells  in  the  three  orthogonal  directions.  This  was 
followed  by  the  affine  distortion  of  the  computational  cell  to 
obtain  the  correct  dimensions  and  density  of  the  computational 
cell;  (b)  a  stochastic  bond-switching  algorithm  was  then 
implemented  (Ref  24)  within  a  Monte  Carlo  computational 
procedure  to  convert  the  crystalline  structure  into  a  disordered 
structure.  Within  this  algorithm,  two  neighboring  Si-0  pairs 
were  randomly  selected  from  the  computational  cell,  and  the 
energy  A E  resulting  from  the  Si-0  bond  switching  was 
computed.  In  the  A E  <  0  case,  the  bond  switching  in  question 
was  accepted  without  any  additional  condition.  On  the  contrary, 
in  the  A E  >  0  case,  a  Boltzmann  probability  factor, 
PB  =  exp(— AE7(3M772))  (where  k  is  the  Boltzmann’s  constant), 
was  first  calculated  and  compared  with  a  random  number  RN 
drawn  from  a  (0,1)  uniform  distribution  function.  The  bond 
switching  in  question  was  then  adopted  only  if  PB  >  RN;  (c)  the 
resulting  structure  was  then  subjected  to  a  carefully  devised  set  of 
NVT  (where  7V(=4608)  is  the  (fixed)  number  of  atoms  within  the 
computational  cell,  V,  the  computational  cell  volume  (also  fixed), 
and  T  is  a  fixed  temperature)  molecular-dynamics  simulations. 
Specifically,  the  NVT  simulations  were  started  at  a  temperature  of 
5300  K  and  carried  out  in  such  a  way  that  the  temperature  was 
controlled  using  the  following  “simulated-annealing”  scheme: 
(i)  a  particle-velocity  scaling  algorithm  was  applied  to  every  time 
step  for  the  first  6000  steps.  This  enforced  strict  control  of  the 
temperature  but  produced  particle  velocities  which  were  incon¬ 
sistent  with  the  target-temperature  Maxwell-Boltzmann  distri¬ 
bution  function;  (ii)  within  the  next  6000  NVT  simulation  steps, 
the  frequency  of  particle-velocity  scaling  was  decreased  to  every 
40  time  steps  while  a  Nose-Hoover  (Ref  25)  temperature  control 
algorithm  (“thermostat”)  was  applied  between  the  particle- 
velocity  scaling  steps.  A  brief  description  of  the  Nose-Hoover 
thermostat  is  provided  in  Appendix  A;  and  (iii)  during  the  final 
8000  steps,  the  temperature  was  controlled  using  only  the  Nose- 
Hoover  thermostat.  This  procedure  ensured  an  efficient  temper¬ 
ature  control  while  yielding  an  equilibrium  state  of  the  material 
(i.e.,  a  particle-velocity  distribution  consistent  with  the  target- 
temperature  Maxwell-Boltzmann  distribution  function). 

Upon  establishing  the  thermodynamic  equilibrium  at 
5300  K,  the  target  temperature  was  reduced  by  500  K,  and 
then  this  procedure  was  re-applied  at  progressively  (by  500  K) 
lower  temperatures,  until  the  final  temperature  of  300  K  was 
reached.  The  total  system  equilibration  procedure  typically 
involved  simulation  times  on  the  order  of  500  ps  resulting  in  an 
average  cooling-rate  of  ~10  K/ps.  An  example  of  a  typical 
molecular-level  topology  within  a  computational  cell  is  dis¬ 
played  in  Fig.  2(a).  A  close-up  of  the  molecular-level  random- 
network  is  displayed  in  Fig.  2(b). 

3.2  Force  Fields 

The  behavior  of  a  material  system  at  the  molecular-level  is 
governed  by  the  appropriate  force  fields  which  describe,  in  an 
approximate  manner,  the  various  interactions  taking  place 
between  the  constituent  particles,  atoms,  ions,  charge  groups, 
etc.  In  other  words,  the  knowledge  of  force  fields  enables 
determination  of  the  potential  energy  of  a  system  in  a  given 
configuration.  In  addition,  gradients  of  the  force-field  functions 
quantify  the  net  forces  experienced  by  the  particles,  the 
information  that  is  needed  in  the  molecular-dynamics 
simulations. 


826— Volume  21(6)  June  2012 


Journal  of  Materials  Engineering  and  Performance 


Fig.  2  (a)  The  computational  unit  cell  for  fused  silica  molecular-le¬ 

vel  simulations  used  in  the  present  study  and  (b)  an  example  of  the 
local  atomic  structure 

In  general,  the  potential  energy  of  a  system  of  interacting 
particles  can  be  expressed  as  a  sum  of  the  valence  (or  bond), 
^valence,  cross-term,  ^cross-term,  and  non-bond,  ^non-bond,  inter- 
action  energies  as 

Aotal  =  ^valence  "E  Across— term  T  ^non— bond  (Eql) 

The  valence  energy  generally  accounts  for  the  contribution 
of  valence  electrons  bonding  and  contains  the  following 
components:  (a)  a  bond  length/stretching  term;  (b)  a  two-bond 
angle  term;  (c)  a  three-bond  dihedral/torsion  angle  term;  (d)  an 
inversion  (or  a  four-atom  out-of-plane  interaction)  term;  and  (e) 
the  so-called  three-atom  Urey-Bradlay  term. 

The  cross-term  interacting  energy  accounts  for  the  cross¬ 
interactions  between  the  aforementioned  valence-energy  com¬ 
ponents  and  includes  terms  like  (a)  stretch-stretch  interactions 
between  two  adjacent  bonds;  (b)  stretch-bend  interactions 
between  a  two-bond  angle  and  one  of  its  bonds;  (c)  bend-bend 
interactions  between  two  valence  angles  associated  with  a 
common  vertex  atom;  (d)  stretch-torsion  interactions  between  a 
dihedral  angle  and  one  of  its  end  bonds;  (e)  stretch-torsion 
interactions  between  a  dihedral  angle  and  its  middle  bond;  (f) 
bend-torsion  interactions  between  a  dihedral  angle  and  one  of 
its  valence  angles;  and  (g)  bend-bend- torsion  interactions 
between  a  dihedral  angle  and  its  two  valence  angles. 

The  non-bond-interaction  term  accounts  for  the  interactions 
between  non-bonded  atoms  and  includes  the  van  der  Waals 
energy  and  the  Coulomb  electrostatic  energy. 

In  the  present  study,  the  so-called  COMPASS  (Condensed- 
phase  Optimized  Molecular  Potentials  for  Atomistic  Simulation 
Studies)  force  field  is  used  (Ref  26,  27).  This  highly  accurate 
force  field  is  of  an  ab-initio  type  since  most  of  its  parameters 
were  determined  by  matching  the  predictions  made  by  the 
ab-initio  quantum  mechanics  calculations  to  the  condensed- 
matter  experimental  data.  A  summary  of  the  COMPASS  force- 
field  functions  can  be  found  in  our  previous  study  (Ref  28). 


3.3  Computational  Method(s) 

All  the  molecular-level  (i.e.,  equilibrium  and  non-equilib¬ 
rium  molecular  dynamics)  calculations  were  carried  out  within 
the  present  study  using  Discover  (Ref  29)  (an  atomic  simulation 
program  from  Accelrys). 

Within  the  equilibrium  molecular-dynamics  method,  the 
system  under  consideration  is  coupled  to  an  (external)  envi¬ 
ronment  (e.g.,  constant-pressure  piston,  constant-temperature 
reservoir,  etc.)  which  ensures  that  the  system  remains  in 
equilibrium  (i.e.,  the  system  is  not  subjected  to  any  thermody¬ 
namic  fluxes).  In  the  present  study,  ATT,  NPT  (P  is  pressure), 
and  NVE  {E  is  the  total  energy)  equilibrium  molecular- 
dynamics  simulations  were  used.  Equilibrium  molecular- 
dynamics  calculations  enable  determination  of  the  (equilibrium) 
thermodynamic  properties  of  a  material  system  through  the  use 
of  time  averages  of  the  state  variables  sampled  along  the 
calculated  system  trajectories. 

Within  non-equilibrium  molecular  dynamics,  the  system  is 
subjected  to  large  perturbations  which  create  a  thermodynamic 
flux  (e.g.,  the  flux  of  energy  and  momentum). 

3.4  Generation  of  Planar  Shock-Waves 

To  generate  a  planar  shock  (or  more  precisely  a  pair  of 
planar  shocks)  within  the  thermodynamically  equilibrated 
material  system  (discussed  earlier),  the  computational  cell 
was  continuously  contracted  along  the  axial  x-direction.  This 
was  accomplished  by  varying/reducing  the  lattice  parameter  a 
as 

a(t)  —  a(t  —  0)  —  2upt  (Eq  2) 

where  t  denotes  time,  up  is  the  so-called  piston  velocity  (or 
equivalently  the  particles  upstream/behind-the-shock  velocity) 
in  the  x-direction.  up  is  varied  over  a  range  between  250  and 
3250  m/s  to  simulate  the  generation  and  propagation  of  shock 
of  various  strengths.  Meanwhile,  computational-cell  trans¬ 
verse  lattice  parameters  b  and  c  are  kept  constant  to  obtain 
planar  (uniaxial-strain)  shock  conditions.  In  this  process,  the 
computational  cell  faces  normal  to  the  shock  propagation- 
direction  behave  very  similarly  to  the  impact-surface  of  a 
plate-like  target  subjected  to  a  so-called  symmetric  “flyer- 
plate”  impact  test  (Ref  30).  Owing  to  the  use  of  the  periodic 
boundary  conditions  in  the  computational  cell  contraction 
direction,  the  procedure  employed  here  generates  a  pair  of 
planar  longitudinal  shocks  which  propagate  at  a  shock  speed 
Us  from  the  computational  cell  faces  (normal  to  the  computa¬ 
tional  cell  contraction  direction)  towards  the  cell  center.  As 
schematically  shown  in  Fig.  3,  these  shocks  leave  behind  a 
“shocked”  material  state  characterized  by  a  higher  material 
density  (as  well  as  by  higher  levels  of  internal  energy,  tem¬ 
perature,  stress,  particle  velocity,  and  entropy). 

3.5  Problem  Formulation 

The  problem  addressed  in  the  present  study  involved 
generation  of  shocks  of  different  strengths  (using  the  afore¬ 
mentioned  computational  cell  parameter  contraction  method), 
determination  of  the  associated  shock-Hugoniot  relations  and 
identification  and  elucidation  of  the  main  molecular-level 
inelastic-deformation/energy  dissipation  processes  taking  place 
at  or  in  the  vicinity  of  the  shock  front.  The  procedure  for  shock- 
wave  generation  was  presented  in  the  previous  section  and  will 
not  be  considered  any  further. 
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Fig.  3  A  schematic  of  the  generation  of  a  pair  of  shocks  in  a 
molecular-level  system  via  the  process  of  computational-cell  parame¬ 
ter  contraction 

As  far  as  the  shock  Hugoniot  determination  is  concerned,  it 
entailed  the  knowledge  of  the  shock-wave  profiles  (and  their 
temporal  evolution)  for  the  axial  stress,  material  density, 
particle  velocity,  internal  energy  and  temperature.  The  latter 
are  obtained  by  lumping  particles/atoms  and  their  (bond  and 
non-bond)  potential  and  kinetic  energy  contribution,  into  fixed- 
width  bins,  in  the  order  of  their  axial  coordinates.  As  will  be 
shown  in  the  next  section,  two  types  of  bins  are  used:  (a)  a 
Lagrangian-type  which  is  fixed  to  the  initial/reference  state  of 
the  computational  cell  and  (b)  a  moving-type  which  is  attached 
to  the  advancing  shock  front. 

Identification  of  the  molecular-level  inelastic-deformation/ 
energy  dissipation  processes  entailed  a  close  examination  of  the 
changes  in  a  material  bond  structure  and  topology  caused  by 
the  passage  of  the  shock. 


4.  Results  and  Discussion 

4.1  Validation  of  the  Computed  Equilibrium  Material  State 

In  this  section,  an  attempt  is  made  to  validate  the  fused  silica 
room  temperature/ambient  pressure  equilibrium  structure 
obtained  through  application  of  the  previously  described 
simulated-annealing  computational  procedure.  In  particular, 
the  material  density,  partial  radial  distribution  functions  for  the 
three  (Si-Si,  Si-O,  and  O-O)  atomic  pairs  and  the  (equilibrium) 
negative  axial-stress  versus  specific  volume  shock-Hugoniot 
relation  are  calculated  and  compared  with  their  experimental 
counterparts. 

To  compute  the  material  density  and  the  partial  radial 
distribution  functions,  NPT  equilibrium  molecular  dynamic 
simulations  were  run  at  the  room  temperature  and  the  ambient 
pressure.  Simulations  were  carried  out  for  over  6000  (0. 1  fs)  time 
steps,  while  maintaining  the  temperature  using  the  Nose-Hoover 
thermostat  and  scaling  the  particle  velocities  for  every  ten  time 
steps.  Pressure  was  controlled  using  a  Berendsen  barostat  (Ref 
3 1).  A  brief  overview  of  this  barostat  is  provided  in  Appendix  B. 
To  obtain  the  equilibrium  Hugoniot  relations,  a  series  of  NVT 


simulations  (described  below)  were  carried  out.  Each  simulation 
was  run  for  over  6000  (0. 1  fs)  time  steps,  while  the  temperature 
was  controlled  using  the  Nose-Hoover  thermostat. 

4.1.1  Density.  The  average  mass-density  of  fused  silica 
computed  from  the  computational  cell  average  volume  and  the 
cell  mass  has  been  found  to  be  1. 5-2.0%  greater  than  the  target 
density  of  2.2  g/cm3.  This  computation/experiment  agreement 
has  been  deemed  to  be  fair. 

4.1.2  Radial  Distribution  Functions.  The  partial  radial- 
distribution  (often  also  referred  to  as  the  partial  pair-correlation) 
function  provides  a  measure  of  the  probability  that,  given  the 
presence  of  an  atom  of  type  a  at  the  origin  of  an  arbitrary 
reference  frame,  there  will  be  an  atom  of  type  /i  within  a 
spherical  shell  of  infinitesimal  thickness  dr  at  a  distance  r  from 
the  reference  atom.  In  amorphous  materials  like  fused  silica,  the 
partial  pair-correlation  functions  are  quite  important  since  they 

(a)  provide  an  insight  into  the  short-range  order  of  the  system; 

(b)  can  be  used  in  the  assessment  of  continuum-level  thermo¬ 
dynamic  material  properties;  and  (c)  provide  a  way  of 
validating  the  molecular-level  calculations  since  these  quanti¬ 
ties  can  be  determined  experimentally  using  X-ray  diffraction. 

The  computed  partial  radial-distribution  functions  are  com¬ 
pared  with  their  counterparts  based  on  the  shell  model 
molecular-dynamics  calculations  (Ref  32).  The  results  of  this 
comparison  are  displayed  in  Fig.  4(a)  and  (b).  Careful  exam¬ 
ination  of  the  present  results  displayed  in  Fig.  4(a)  shows  that 
they  are  qualitatively  similar  to  the  ones  reported  in  Ref  32.  As 
far  as  the  quantitative  agreement  between  the  two  sets  of  results 
is  concerned,  it  could  be  characterized  as  being  fair. 

4.1.3  Negative  Axial-Stress  versus  Specific  Volume  Equi¬ 
librium  Shock  Hugoniot.  Finally,  the  method  of  Erpenbeck 
(Ref  33)  is  employed  to  compute  the  “equilibrium”  negative 
axial-stress  versus  specific  volume  Hugoniot  relations  for  fused 
silica.  Toward  that  end,  a  Rankine-Hugoniot  function  is  first 
constructed  as 


H(T,  V)  =E(T,  V)-Eo(To,V0) 

+  \(~tn(T,  V)  +  -tn,0(T,  V))(V  -  V0) 

(Eq  3) 

where  tn  denotes  the  axial  stress,  and  the  subscript  “0” 
denotes  the  initial  (unshocked)  material  state.  In  order  to  com¬ 
pute  a  single  point  of  the  negative  axial-stress  versus  specific 
volume  Hugoniot,  a  series  of  NVT  equilibrium  molecular- 
dynamics  simulations  was  carried  out.  In  each  of  these  series 
of  simulations,  the  computational-cell  volume  was  kept  con¬ 
stant  while  the  temperature  was  varied  between  different  simu¬ 
lations  within  the  same  series.  Each  simulation  within  a  given 
(fixed  volume)  series  yielded  the  associated  negative  axial 
stress  and  energy  density  values.  These  were  next  used  to 
evaluate  the  Rankine-Hugoniot  function,  Eq  3  at  each  simula¬ 
tion  temperature.  Then  the  value  of  the  temperature  at  which 
the  Rankine-Hugoniot  function  is  equal  to  zero  is  determined 
via  interpolation.  Finally,  an  NVT  run  is  carried  out  at  this 
temperature  to  determine  the  corresponding  negative  axial- 
stress  (and  energy  density).  This  procedure  is  then  repeated 
for  a  set  of  (fixed)  computational  cell  volumes,  to  construct 
the  negative  axial-stress  versus  specific  volume  Hugoniot 
(complete)  relation.  The  —tn  versus  v/v0  Hugoniot,  obtained 
using  the  procedure  described  above,  is  depicted  in  Fig.  5.  For 
comparison,  the  corresponding  Hugoniot  curve  obtained 
experimentally  in  Ref  34  is  also  shown  in  this  figure. 


828— Volume  21(6)  June  2012 


Journal  of  Materials  Engineering  and  Performance 


0  20  40  60  SO 
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Fig.  4  A  comparison  between  the  fused  silica  Si-Si,  Si-O,  0-0  par¬ 
tial  radial-distribution  functions  computed  in  (a)  the  initial  equilib¬ 
rium  state,  as  computed  in  the  present  study;  (b)  the  equilibrium 
initial  state  as  reported  in  Tilocca  et  al.  (Ref  32);  and  (c)  the 
as-shocked  material  state,  as  computed  in  the  present  study 

A  comparison  between  the  computed  and  the  experimental 
(Ref  3)  Hugoniot  curves  displayed  in  Fig.  5  shows  that  the 
agreement  is  relatively  good  at  small  compression  levels  (the 
weak  shock  regime)  and  it  worsens  as  the  shock  strength  is 
increased  (i.e.,  v/v0  is  decreased).  Specifically,  at  the  lowest 


Specific  Volume  over  Zero-Pressure  Specific  Volume,  v/v0 


Fig.  5  A  comparison  between  the  computed  equilibrium,  the  com¬ 
puted  non-equilibrium  and  the  experimentally  measured  (Ref  3)  neg¬ 
ative  axial-stress  versus  normalized  specific  volume  Hugoniot  curves 
in  fused  silica 


values  of  v/v0,  the  negative  axial  stress  is  under-predicted  by 
35-40%. 

It  should  be  noted  that  the  Hugoniot  relations  obtained  in 
this  portion  of  the  work  are  referred  to  as  being  of  the 
“equilibrium”  character.  This  was  done  to  emphasize  the  fact 
that  these  relations  were  obtained  using  equilibrium  molecular- 
dynamics  simulation  results.  In  other  words,  the  determination 
of  the  Hugoniot  relations  was  carried  out  under  the  assumption 
that  shock  loading  transitions  the  initial  ambient-pressure 
equilibrium  material  state  into  another,  high-pressure,  equilib¬ 
rium  material  state,  (i.e.,  the  dissipative  portion  of  the  internal 
energy  change  accompanying  the  passage  of  a  shock  is  purely 
thermal).  However,  shock  loading  is  not  only  an  intrinsically 
irreversible  non-equilibrium  energy-dissipating  process,  but 
may  also,  in  the  case  of  solid  materials,  be  accompanied  by 
various  microstructure-altering  processes  (e.g.,  extensive 
inelastic  deformations,  damage,  phase  transformations,  etc.). 
The  operation  of  such  processes  would  cause  the  as-shocked 
material  state  not  to  lie  any  longer  on  the  equilibrium  equation- 
of-state  surface  and,  hence,  one  may  question  the  validity  of  the 
equilibrium  Hugoniot  relations  (particularly  in  the  regime 
associated  with  intermediate  and  strong-shocks  (the  regime  in 
which  the  equilibrium  Hugoniot  relations  are  found  not  to  agree 
well  with  their  experimental  counterparts,  Fig.  5)).  To  test  the 
validity  of  these  relations,  non-equilibrium  molecular-dynamics 
simulations  are  used  to  generate  and  drive  planar  shocks  in 
fused  silica.  This  provided  an  alternative,  more-direct  way  for 
determination  of  the  shock-Hugoniot  relations.  Detailed  results 
related  to  this  portion  of  the  study  are  presented  in  the 
remainder  of  the  article. 

4.2  Shock-Wave  Molecular-Level  Analysis 

As  explained  earlier,  the  computational  cell  contraction 
procedure  yielded  a  pair  of  converging  shocks  per  computa¬ 
tional  cell  (due  to  the  use  of  the  periodic  boundary  conditions). 
Since  the  two  shocks  were  structurally  identical  and  the 
collision  of  the  two  converging  shocks  is  beyond  the  scope  of 
the  present  study,  the  analysis  presented  in  this  and  the 
subsequent  sections  focus  on  the  behavior  of  the  right 
propagating  shock  only. 
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4.2.1  Identification  of  the  Shock-Wave  Front.  Fig¬ 
ure  6(a)  and  (b)  show  examples  of  the  typical  results  obtained 
in  this  portion  of  the  study.  These  results  pertain  to  material 
molecular-level  microstructure/topology  evolution  accompany¬ 
ing  the  propagation  of  a  pair  of  approaching  shocks.  Arrows  are 
used  in  Fig.  6(a)  and  (b)  to  indicate  an  approximate  location  of 
the  mid-plane  of  the  shock  front.  The  shock  front  was 
determined  as  the  region  within  which  a  rapid  change  in  the 
material  density  is  detected.  A  careful  examination  of  the 
results  displayed  in  these  figures  shows  that  the  shock  remains 
fairly  planar  during  its  motion.  Clearly,  this  finding  was 
influenced  by  the  fairly  limited  computational-cell  size  in  the 
lateral  direction  as  well  as  by  the  imposition  of  the  periodic 
boundary  conditions. 

4.2.2  Shock- Wave  Front  Analysis.  To  obtain  a  more 
quantitative  description  of  the  shock-wave  form,  a  set  of 
previously  described  parallel  Lagrangian  bins  stacked  in  the 
shockwave  propagation  direction  was  utilized.  As  explained 
earlier,  these  bins  are  fixed  to  the  reference  frame  of  the 
material  and,  hence,  contain  the  same  set  of  atoms  throughout 
the  simulation.  This  enables  the  monitoring  of  the  average 
properties  of  fixed  sets  of  atoms  to  reveal  the  structure  of  the 
advancing  shock. 

Figure  7(a)  and  (b)  display  the  results  obtained  through  the 
use  of  the  aforementioned  method  in  which  x-component  (i.e., 
shock  propagation  direction)  particle  velocities  at  different 


Fig.  6  Temporal  evolution  of  the  molecular  level  material  micro¬ 
structure  accompanying  generation  and  propagation  of  a  pair  of  pla¬ 
nar  shocks  in  fused  silica.  Arrows  denote  the  position  of  the  shock 
front  mid-plane 


simulation  (i.e.,  post-shock-wave  generation)  times  are  plotted 
against  the  Lagrangian-bin  mid-plane  x-location.  The  results 
displayed  in  Fig.  7(b)  are  obtained  under  identical  simulation 
conditions  except  for  the  rate  of  axial  contraction  of  the 
computational  cell  (3000  m/s  in  Fig.  7(a)  and  2400  m/s  in 
Fig.  7(b)). 

A  brief  examination  of  the  results  displayed  in  Fig.  7(a)  and 
(b)  reveals  that 

(a)  two  shocks  are  generated  (only  the  right-propagating 
shock  is  shown,  though)  at  the  computational  cell  faces 
normal  to  the  x-direction.  These  shocks  subsequently 
propagate  towards  the  computational-cell  center; 

(b)  following  a  brief  transient  period,  the  (right-propagating) 
shock  appears  to  become  steady  (i.e.,  within  a  reference 
frame  which  is  attached  to,  and  moves  with,  the  shock- 
wave  front,  the  wave  front  is  nearly  time-invariant); 

(c)  the  particle  velocity  and  the  shock  speed  increase  while 
shock-wave  width  decreases  with  an  increase  in  the  com¬ 
putational-cell  contraction  rate.  It  should  be  noted  that  the 
curves  displayed  in  Fig.  7(a)  and  (b)  are  labeled  using  the 
associated  post-shock  generation  times  in  femtoseconds. 


(a)  Distance  in  x-direction,  nm 


(b)  Distance  in  x-direction,  nm 


Fig.  7  Temporal  evolution  of  the  particle  velocity  associated  with 
the  propagation  of  a  right-propagating  planar  shock  in  fused  silica 
under  the  imposed  computational  cell  contraction  rate  of  (a)  3000  ml 
s  and  (b)  2400  m/s.  Curve  labels  denote  the  associated  post  shock- 
generation  time  in  femtoseconds 
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It  should  also  be  noted  that  no  (artificial  viscous-damping 
based)  thermostat  was  used  in  the  present  non-equilibrium 
molecular-dynamics  simulations.  Thus,  the  observed  steady- 
shock-wave  profile  is  a  natural  consequence  of  a  balance 
between  the  continuous  supply  of  energy  to  the  system 
(through  the  continuous  computational  cell  axial  contraction) 
and  the  attendant  bond/microstructure-altering  atomic  motions 
in  the  continuously  enlarged  upstream  material  region  swept  by 
the  shock.  Further  examination  of  the  results  displayed  in 
Fig.  7(a)  and  (b)  reveals  that  in  the  shockstrength  range 
examined,  the  observed  shock  width  is  quite  small  (ca.  2-3  nm) 
and  that  it  decreases  with  an  increase  in  shock  strength.  This 
finding  suggests  that  the  dominant  energy  dissipation  processes 
captured  by  the  molecular-level  simulations  of  shock  genera¬ 
tion/propagation  within  fused  silica  are  quite  weak.  As  will  be 
shown  later,  these  processes  involve  changes  in  the  fused  silica 
atomic  coordination  and  in  the  size  of  the  smallest  Si-0  rings. 
In  other  words,  the  known  viscous  dissipation  processes  with 
characteristic  times  at  the  microsecond  time  scale  are  not  active 
in  (captured  by)  the  present  simulations.  Based  on  the 
fundamental  analysis  of  propagation  of  the  structured  longitu¬ 
dinal  planar  shocks  (Ref  35),  this  fact,  in  general,  only  affects 
the  width  of  the  shock  profile  but  not  the  values  of  the  as- 
shocked  material  state  variables. 

4.2.3  Non-Equilibrium  Shock  Hugoniot  Relations.  The 
results  like  the  ones  displayed  in  Fig.  7(a)  and  (b)  were  utilized 
to  determine  a  functional  relationship  between  the  (Lagrangian) 
shock  speed,  Us,  and  the  particle  velocity,  up.  The  outcome  of 
this  procedure  is  displayed  in  Fig.  8.  The  Us  versus  up  relation 
like  the  one  displayed  in  Fig.  8  is  generally  considered  as  one 
of  the  shock  Hugoniot  relations.  Since,  as  discussed  earlier,  the 
Hugoniot  is  geometrically  a  hyper-line  in  a  multi-dimensional 
stress/pressure,  energy  density,  mass  density  (or  specific 
volume),  temperature,  particle  velocity,  shock-speed  etc.  space. 
The  Us  versus  up  (Fig.  8)  and  the  negative  axial-stress  versus 
specific  volume  (Fig.  5)  relations  are  simply  the  particular 
planar  (two-dimensional)  projections  of  the  Hugoniot.  The  Us 
versus  up  Hugoniot  data  displayed  in  Fig.  8  are  fitted  to  a  linear 
relationship  in  the  form  Us  =  cQ  +  sup ,  and  the  following  values 
of  the  two  coefficients  are  obtained  as  c0  =  2560  m/s  and 
s  =  1.182. 


Particle  Velocity  xdot  (m/s) 

Fig.  8  Shock-speed  vs.  particle-velocity  Hugoniot  data  (based  on 
the  present  results  like  the  ones  displayed  in  Fig.  7(a)  and  (b))  and 
the  corresponding  linear  fit  to  the  data 


The  other  commonly  used  planar  shock  Hugoniot  relations 
include  negative  axial-stress  versus  mass-density,  p;  (mass- 
based)  internal  energy  density,  E ,  versus  p  (or  v);  —tn  versus 
mp,  and  temperature  T  versus  p  (or  v).  Following  our  prior 
studies  (Ref  9,  21),  an  attempt  was  made  to  determine  these 
relations  using  two  distinct  methods: 

(a)  A  method  based  on  the  Lagrangian  form  of  the  three 


jump  equations  defined  as 

P  =  p  +(US  —  i/p) 

(Eq  4) 

hi  +  P  ~  *11  +  P+(^s  —  Mp) 

(Eq  5) 

e~  +  +  0.5 Ul  =  e+  +  +  0.5  (Us  -  up)2 

(Eq  6) 

The  above  equations  relate  the  known  downstream  material 
states  (denoted  by  a  superscript  “  — ”)  and  the  unknown 
upstream  material  states  (denoted  by  a  superscript  “+”) 
associated  with  the  shock  of  a  given  strength  (as  quantified 
by  the  shock  speed  or  the  downstream-to-up stream  particle 
velocity  jump).  These  equations  are  next  combined  with  the 
previously  determined  Us  versus  up  relation  and  the  prescribed 
(shock-strength  defining  quantity)  Us  or  up  to  solve  for  the 
unknown  upstream/as-shocked  material  states.  It  should  be 
noted  that  this  method  enables  determination  of  only  material 
mechanical  state  variables  (tn,  E ,  v(=l/p),  Us  and  up).  To  obtain 
temperature-based  Hugoniot  relations,  the  procedure  described 
recently  in  Ref  36  was  employed.  This  procedure  involves  the 
following  steps:  (i)  the  (differential)  statement  of  the  combined 
first  and  second  laws  of  thermodynamics  is  equated  to  the 
differential  form  of  the  Rankine-Hugoniot  equation;  (ii)  the 
term  representing  the  rate  of  change  of  mass-based  entropy 
density  with  the  specific  volume  appearing  in  the  resulting 
equation  is  eliminated  using  the  expression  for  the  total 
derivative  of  temperature  and  the  Maxwell  equation  based 
thermodynamic  relations;  and  (iii)  the  resulting  differential 
equation  which  now  contains  additional  material  parameters 
(the  constant-volume  specific  heat  and  the  Gruneisen  Gamma) 
and  the  stress/pressure  Hugoniot  relation  is  integrated  numer¬ 
ically  between  the  initial  and  the  final  material  states  to  obtain 
the  temperature  corresponding  to  the  given  shock  strength;  and 
(b)  Thermodynamic  averages  of  the  as-shocked  material  in  the 
shock  are  used  to  calculate  the  temperature  (as  well  as  the 
remaining  material  state  variables  like  stress/pressure,  specific- 
volume/density,  mass-based  energy  density,  etc.).  The  neces¬ 
sary  data  were  collected  using  non-Lagrangian  bins.  These  bins 
were  attached  to  (and  moved  with)  the  steady-shock  front  and 
hence,  in  contrast  to  the  Lagrangian  bins,  collected  the 
information  about  the  atoms  (temporarily)  residing  in  a  given 
segment  of  the  steady-shock  profile  rather  than  the  information 
about  a  fixed  set  of  atoms.  Also,  since  the  main  role  of  these 
bins  was  to  collect  the  data  in  the  upstream  material,  they  were 
placed  only  in  the  regions  behind  the  shock.  This  procedure 
was  found  to  be  quite  work  intensive  and  to  yield  widely 
scattered  results.  To  overcome  these  shortcomings,  an  attempt 
was  made  to  carry  out  NVE  (E  is  the  total  energy  of  the 
computational  cell)  molecular-dynamics  simulations  (over  an 
extended  time  period)  on  the  material  extracted  from  the 
upstream  region  and  placed  into  its  own  periodic  cell. 
Unfortunately,  this  attempt  was  not  very  successful  since 
prolonged  molecular  level  simulations  of  the  “shocked  mate¬ 
rial”  were  associated  with  extensive  material  relaxation  so  that 
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the  material  molecular  structure  did  not  any  longer  correspond 
to  the  as-shocked  material  state.  Consequently,  the  shock- 
Hugoniot  results  obtained  using  this  procedure  are  deemed 
unreliable  and  are  not  shown  and  not  used  in  the  remainder  of 
this  article. 

The  results  of  the  method  (a),  based  on  the  use  of  the  shock- 
jump  equations,  are  depicted  in  Fig.  9(a)-(d).  These  four  figures 
show  the  —tn  versus  p,  E  versus  p,  —tn  versus  uv ,  and  T 
versus  p  Hugoniot  relations,  respectively.  In  addition,  the  non¬ 
equilibrium  —tn  versus  v/v0  Hugoniot  relation  is  depicted  in 
Fig.  5  and  compared  with  its  equilibrium  computational  and 
experimental  counterparts.  It  should  be  noted  that  only  the  non¬ 
equilibrium  molecular-dynamics  results  associated  with  the 
strong-  and  moderate-shock  regimes  are  displayed  in  Fig.  5. 
The  reason  for  this  is  that,  under  weak-shock  conditions,  the 
results  did  not  show  the  formation  of  a  well-defined  shock.  This 
finding  is  consistent  with  the  experimental  Hugoniot  curve  in 
Fig.  5  which,  in  the  weak-shock  regime,  is  concave  downward, 
suggesting  an  abnormal  (i.e.,  a  non-shock-supporting)  behavior 
of  fused  silica  under  compressive  loading.  In  other  words,  in 
the  weak-shock  regime,  compressive  loading  tends  to  produce 
spreading  finite-amplitude  waves  rather  than  shocks.  Examina¬ 
tion  of  Fig.  5  reveals  that  the  non-equilibrium  computed 
Hugoniot  data  are  in  substantially  better  agreement  with  their 


experimental  counterparts  than  are  the  equilibrium  computed 
results. 

The  Hugoniot  relations  displayed  in  Fig.  5  and  Fig.  9(a)-(d) 
are  typically  used  within  a  continuum-level  computational 
analysis  of  shock-wave  generation/  propagation  in  two  ways: 

(a)  They  are  directly  used  in  the  analysis  of  planar  longitu¬ 
dinal  shock-wave  propagation  under  uniaxial  strain  con¬ 
ditions  (Ref  37)  and 

(b)  Alternatively,  they  can  be  used  to  derive  a  continuum-le¬ 
vel  material  model  which  is  consistent  with  the  material 
mechanical  response  under  high-rate,  large-strain,  high- 
pressure  conditions.  Such  a  model  is  subsequently  used 
in  general  three-dimensional,  nonlinear  dynamics  com¬ 
putational  analyses  (Ref  38,  39). 

4.2.4  Shock-Induced  Material  Microstructure-Evolu- 
tion.  The  results  presented  and  discussed  in  the  previous 
sections  clearly  revealed  the  formation  and  propagation  of 
planar  shocks  in  fused  silica  and  enabled  formulation  of  the 
appropriate  shock-Hugoniot  relations.  In  the  present  section,  a 
more  detailed  investigation  of  the  molecular-level  material 
micro  structure,  in  the  wake  of  a  propagating  planar  shock,  is 
carried  out. 


Fig.  9  (a)  Negative  axial-stress  vs.  mass-density;  (b)  internal  energy  density  vs.  mass-density,  (c)  negative  axial  stress  vs.  particle  velocity,  and 

(d)  temperature  vs.  mass-density  Hugoniot  relations  for  fused  silica  obtained  using  non-equilibrium  molecular-dynamics  calculations 
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At  the  onset,  it  should  be  noted  that,  the  absence  of  a  crystal 
structure  in  fused  silica  makes  the  analysis  of  material 
micro  structure  and  its  evolution  a  formidable  task.  That  is, 
the  presence  of  a  crystal  lattice  (as  is  the  case  in  single-crystal 
solids)  studied  earlier  (Ref  40,  41)  makes  the  analysis  of  shock¬ 
loading  induced  changes  and  the  associated  deformation 
processes  in  the  material  micro  structure  more  tractable.  Fused 
silica,  on  the  other  hand,  is  an  amorphous  material  in  both  its 
initial  and  its  as-shocked  state.  To  quantify  the  extent  of  shock- 
induced  changes  in  the  microstructure  of  fused  silica,  the 
following  microstructural  parameters  were  tracked:  (a)  the 
random  network  X  parameter,  defined  earlier.  The  R  parameter 
(=2.0)  was  not  monitored  since  its  value  is  controlled  by 
glass-chemistry  and  not  by  the  material  microstructure.  The 
^parameter  was  not  tracked  either  since  in  fused  silica  the  X  +  Y 
sum  is  equal  to  4.0;  (b)  the  size-distribution  of  the  smallest  Si-0 
rings;  (c)  the  Sz-atom  average  coordination  number  (i.e., 
bonding  structure);  and  (d)  the  radial  distribution  function. 

Changes  in  X  Parameter.  Shock-induced  changes  in  the 
random-network  X  parameter  are  found  to  be  the  result  of  a 
competition  between  Si-0  bond  breaking  processes  (increase 
the  average  number  of  non-bridging  oxygens  per  un-shocked 
material  polyhedron)  and  the  ones  which  lead  to  an  increase  in 
Si-atom  coordination  number  (decrease  the  A-parameter  value). 
The  results  obtained  in  the  present  study  show  that  the  instances 
of  Si-atom  coordination  number  increase  greatly  out-number 
the  ones  associated  with  Si-0  bond-breaking.  Consequently,  the 
value  of  the  X  parameter  has  been  found  to  decrease  by  5-10% 
(with  the  largest  percent  changes  observed  under  the  strongest- 
shock  loading  conditions). 

Size  Distribution  of  the  Smallest  Si-O  Rings.  It  is  customary  to 
monitor  changes  in  the  size  distribution  of  the  smallest  Si-0 
rings  and  use  the  same  as  an  evidence  for  the  attendant  shock- 
induced  microstructural  changes.  However,  one  should  exercise 
caution  at  this  juncture.  It  should  be  first  noted  that  the 
crystobalite  allotropic  modifications  of  Si02  contain  only  the 
so-called  six-membered  rings  (i.e.,  rings  containing  six  silicon 
and  six  oxygen  atoms).  In  quartz,  on  the  other  hand,  one  finds 
both  six  member  and  eight  member  Si-0  rings.  In  both  cases, 
the  total  number  of  the  smallest  rings  associated  with  each  Si 
atom  is  six,  and  silicon  atoms  are  four-fold  coordinated. 
Intuitively,  one  would  expect  that  a  larger  fraction  of  smaller 
rings  in  a  given  allotropic  modification  would  yield  a  material 
with  closer  packing  and,  thus,  a  higher  density.  However,  this  is 
not  the  case  considered  for  this  study  since  the  density  of 
crystobalite  is  2.17  g/cm3  while  that  for  quartz  is  2.64  g/cm3. 
The  structure  of  the  six-membered  rings  in  crystobalite  and 
quartz  are  displayed  in  Fig.  10(a)  and  (b).  Examination  of 
Fig.  10(a)  and  (b)  reveals  that  the  atomic  structures  of  the  six- 
membered  rings  in  the  two  Si02  allotropes  are  quite  different 
and  that  these  rings  in  quartz  occupy  less  space  than  the  ones  in 
crystobalite.  Furthermore,  in  stishovite  in  which  silicon  atoms 
have  six-fold  coordination,  one  finds  the  two-  and  three- 
membered  Si-0  rings  and  the  density  is  4.29  g/cm3.  This  brief 
analysis  suggests  that  there  is  no  simple  correlation  between  the 
size  of  the  smallest  Si-0  rings  (as  measured  by  the  number  of  Si 
and  O  atoms  in  the  rings)  and  the  fused  silica  material  density 
or  its  microstructure.  Instead,  the  information  regarding  the 
size,  topology,  and  distribution  of  the  smallest  rings  should  be 
combined  with  the  silicon  atom  coordination  data. 


Fig.  10  Two  views  of  the  atomic  structure  of  six-membered  Si-0 
rings  in  (a)  crystobalite  and  (b)  quartz 

To  quantify  the  size-distribution  of  the  smallest  Si-0  rings  in 
both  the  initial  and  the  as-shocked  fiised-silica  states,  a 
computational  method  was  developed  as  part  of  the  present 
study.  This  method  solves  the  class  of  so-called  shortest-path 
problems  and  is  a  simple  modification  of  the  Dijkstra’s 
algorithm  (Ref  42).  The  main  modifications  in  this  algorithm 
are  associated  with  the  fact  that,  in  the  present  case,  the  starting 
point  and  the  destination  point  of  the  path  are  identical.  The 
smallest-ring  size-distribution  results  for  the  initial  and  two 
representative  as-shocked  material  states  are  displayed  in 
Fig.  11(a)  and  (b),  respectively. 

Since,  the  as-shocked  results  obtained  under  different 
loading  conditions,  Fig.  11(b),  are  mutually  similar  and  sub¬ 
stantially  different  than  the  ones  associated  with  the  initial  state, 
Fig.  11(a),  one  can  conclude  that  the  material  undergoes 
detectible  microstructural  changes  during  shock  loading.  This 
finding  was  further  supported  by  the  observations  that,  while 
the  six-membered  Si-0  rings  in  the  fiised-silica  initial  state 
were  similar  to  those  found  in  crystobalite,  Fig.  10(a),  the  six- 
membered  Si-0  rings  observed  in  the  as-shocked  material  state 
resembled  more  to  the  ones  found  in  quartz,  Fig.  10(b).  Also, 
no  five-membered  rings  were  found  in  the  initial  states  while 
the  same  were  observed  in  the  as-shocked  state.  An  example  of 
the  five-membered  Si-0  ring  found  in  fused-silica  as-shocked 
state  is  displayed  in  Fig.  12(a).  For  clarity,  the  silicon  and 
oxygen  atoms  associated  with  the  five-membered  ring  are 
colored  white  and  green,  respectively,  while  the  remaining  Si 
and  O  atoms  are  colored  yellow  and  red,  respectively. 

Si- Atom  Average  Coordination  Number.  Examination  of  the 
as-shocked  material  micro  structure  of  fused  silica  revealed  that 
the  Si-atom  average  coordination  number  increases  from  its 
initial  value  of  4.0-4.06.  An  example  of  the  shock-loading- 
induced  increase  in  Si-atom  coordination  number  is  depicted  in 
Fig.  12(b).  In  this  figure,  a  five-fold-coordinated  Si-atom 
(white)  is  displayed  along  with  the  associated  five  oxygen 
atoms  (green).  The  remaining  Si  and  O  atoms  are  displayed  as 
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(a)  Smallest  Si-O  Ring-size 


Fig.  11  Size  distribution  function  for  the  smallest  Si-0  rings  in  the 
fused  silica:  (a)  initial  state  and  (b)  two  as-shocked  states 


yellow  and  red,  respectively.  The  forgoing  finding  regarding  an 
increase  in  the  Si-atom  coordination  is  an  indication  of  the 
transition  of  fused  silica  to  a  “stishovite-like”  state,  an 
energetically  favored  state  at  high  shock-induced  stresses.  This 
finding  is  consistent  with  the  experimental  observation  of 
Alexander  et  al.  (Ref  43)  who  reported  the  formation  of 
stishovite-like  domains  in  various  silica-based  glasses  after 
being  subjected  to  shock  loading  in  excess  of  30  GPa. 

Radial  Distribution  Function.  The  last  microstructure-charac- 
terizing  quantity  covered  in  the  present  analysis  is  the  radial 
distribution  function.  In  Fig.  4(c),  the  Si-Si,  Si-O,  and  0-0 
partial  radial  distribution  functions  are  displayed  for  a  typical 
as-shocked  material  state.  Comparison  of  the  results  displayed 
in  Fig.  4(a)  and  (c)  provides  additional  evidence  that  shock 
loading  gives  rise  to  the  changes  in  fused-silica  material 
microstructure. 

4.3  Final  Remarks 

Previous  non-equilibrium  molecular-dynamics  simulations 
of  shock-wave  propagation  in  single-crystalline  solids  (Ref  40, 
41)  established  that  energy  dissipation  accompanying  shock 
propagation  is  the  result  of  inelastic  deformation  (permanent 


Fig.  12  Examples  of  (a)  a  five-membered  Si-0  ring  and  (b)  five¬ 
fold  coordinated  Si-atom  in  as-shocked  fused-silica  material  state. 
The  silicon  atom  in  question  is  colored  white  (the  remaining  Si  atoms 
are  shown  as  yellow),  and  the  associated  O-atoms  are  shown  as  green 
(the  remaining  O-atoms  are  shown  as  red)  (Color  figure  online) 


slippage  of  crystal  planes  and  the  formation  of  crystal  defects) 
and  phase-transformation  processes  and  not  a  result  of 
(velocity-dependent)  viscous  dissipation  (as  is  the  case  for 
shocks  in  fluids).  The  present  computational  results  obtained 
in  fused  silica  suggest  that  this  may  be  also  true  in  the  case 
of  amorphous  solids.  This  finding  provides  a  possible 
explanation  for  the  fact  that  the  axial-stress  versus  specific 
volume  equilibrium  and  non-equilibrium  Hugoniots  differ 
significantly  under  moderate  and  strong  shock  regimes  and 
the  non-equilibrium  Hugoniot  is  in  better  agreement  with  the 
experimental  one.  In  other  words,  the  equilibrium  Hugoniot 
was  derived  under  the  assumption  that  the  material  micro¬ 
structure  does  not  change  under  shock  loading.  When  such 
changes  take  place,  the  Hugoniot  states  no-longer  lie  on  the 
initial  equation-of-state  surface  but  rather  on  a  surface(s) 
associated  with  the  transformed  material  state.  Since  the 
energy  levels  associated  with  the  transformed-material  equa- 
tion-of-state  surface  are  generally  higher  when  compared  to 
their  initial-material  counterparts,  the  corresponding  axial- 
stress  values  are  expected  to  be  higher.  This  was  indeed 
confirmed  by  the  non-equilibrium  molecular-dynamics  results 
displayed  in  Fig.  5. 

In  the  present  study,  molecular-level  computational  methods 
and  tools  are  employed  to  investigate  shock-wave  generation 
and  propagation  phenomena  in  fused  silica.  However,  it  is  not 
the  objective  here  to  suggest  that  shock- wave  molecular-level 
simulations  should  replace  the  corresponding  hydrodynamic/ 
continuum-level  computations.  In  fact,  the  latter  are  generally 
preferred  in  the  studies  dealing  with  real-world  technical 
systems  (e.g.,  blast  impact  of  vehicle  windshields).  On  the  other 
hand,  molecular-level  simulations  appear  to  be  highly  benefi- 
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cial  in  determining  the  shock-Hugoniot  relations  and  in 
identifying  key  microstructure-altering  and  deformation  pro¬ 
cesses  accompanying  shock  propagation.  As  pointed  out 
earlier,  the  knowledge  of  the  shock  Hugoniot  relations  and 
the  associated  phase-transformation/deformation  processes  can 
be  used  to  derive  physically  based  continuum-level  material 
models  which  are  used  in  the  aforementioned  hydrodynamic 
analyses. 


5.  Summary  and  Conclusions 

Based  on  the  results  obtained  in  the  present  study,  the 
following  summary  remarks  and  main  conclusions  can  be 
drawn: 


Appendix  A:  Nose-Hoover  Thermostat  (Ref  25) 

The  simulated-annealing  computational  procedure  described 
in  section  3.1  employed  NVT  molecular-dynamics  simulations 
in  which  temperature  was  controlled  using  the  Nose-Hoover 
thermostat  (Ref  25). 

In  this  section,  a  brief  description  is  provided  of  this 
temperature-control  algorithm.  Within  this  algorithm,  the  real 
(constant-temperature)  molecular-level  material  system  charac¬ 
terized  by  the  particle  coordinates  q\  and  linear  momenta  p\ 
(i  =  1,  2,. .  .,3  N),  as  well  as  by  a  fixed  temperature  T is  replaced 
with  an  extended  (constant-energy)  system  containing  modified 
degrees  of  freedom  qt  —  q\,  pt  —  j  and  two  additional  degrees 
of  freedom  s  (a  dimensionless  coordinate)  and  ps  (the  conjugate 
linear  momentum).  The  extended  system  Hamiltonian  (the  total 
energy)  is  then  defined  as 


1 .  Equilibrium  and  non-equilibrium  molecular-dynamics 
computational  analyses  are  employed  to  study  the  propaga¬ 
tion  of  planar  longitudinal  shocks  in  fused  silica,  a  material 
commonly  used  in  transparent  armor  applications. 

2.  To  create  molecular-level  shocks,  a  computational  proce¬ 
dure  based  on  continuously  contracting  computational 
cell  is  used.  This  procedure  does  not  require  the  use  of  a 
viscous-dissipation-based  thermostat  to  form  a  steady 
shock  wave. 

3.  By  properly  post-processing  the  computational  results 
pertaining  to  the  atomic  positions,  velocities,  and  interac¬ 
tion  forces,  appropriate  shock  Hugoniot  relations  are  de¬ 
rived.  These  relations  define  the  locus  of  stress,  energy 
density,  mass  density,  temperature,  and  particle  velocity 
in  the  as-shocked  states  of  fused  silica. 

4.  Detailed  examination  of  the  molecular-level  microstruc¬ 
ture  in  the  as-shocked  state  of  fused  silica  is  carried  out 
to  identify  the  nature  of  the  microstructure-altering  and 
deformation  processes  responsible  for  energy  dissipation 
accompanying  shock-wave  propagation.  This  exanimation 
identified  significant  and  shock-strength-dependent 
changes  in  the  silicon  atom  coordination  on  the  structure 
of  the  Si-0  random-network  connectivity.  In  particular,  a 
substantially  larger  number  of  five-fold  coordinated  Si 
atoms  has  been  found  in  the  as-shocked  state  relative  to 
that  found  in  the  pre-shocked  fused-silica  initial  state. 
Also,  smaller  Si-0  rings,  not  found  in  the  initial  state, 
have  been  observed  in  the  fused  silica  with  as-shocked 
state.  It  has  been  subsequently  shown  that  these  micro- 
structural  changes  give  rise  to  a  deviation  of  the  equa- 
tion-of-state  surface  for  the  as-shocked  material 
compared  with  the  initial  material.  Consequently,  higher 
stress-levels  are  attained  in  the  as-shocked  material  than 
the  ones  which  are  predicted  by  an  analysis  based  on  the 
equilibrium  equation-of-state. 
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H  =  Y^P^/ZmiS2  +  q>{qi)  +  p2j2Q  +  gkT\ns  (Eq  A.l) 

i 

while  the  associated  equations  of  motion  are  defined  as 


dqt  dH  pi  2 

dt  dpi  Mi 


(Eq  A.2) 


dpi  dH  d(p 

dt  dqt  dqi 


(Eq  A.3) 


ds  dH  ps 
dt  dps  Q 


(Eq  A.4) 


dps  dH_  (^jP[/2misYzSPl  (Eq  A  5) 

dt  ds  s 

where  Q  is  a  mass  term,  associated  with  the  added  degrees 
of  freedom,  cp(^)  is  the  real-system  force-field  potential, 
g  (=3  N),  and  k  is  the  Boltzmann’s  constant.  The  tempera¬ 
ture  control  is  then  carried  out  by  integrating  the  extended 
system  in  time  and  projecting  the  (computed)  extended 
(NVE  canonical-ensemble  related)  system  trajectory  data 
onto  the  real  ( NVT  microcanonical  ensemble  related) 
system. 


Appendix  B:  Berendsen  Barostat  (Ref  31) 

Within  this  barostat,  pressure  is  controlled  by  scaling  the 
particle  coordinates  and  the  computational-cell  edge  lengths 
while  maintaining  the  shape  of  the  computational  cell.  This,  in 
a  way,  is  similar  to  coupling  the  system  to  a  constant-pressure 
environment.  The  strength  of  the  coupling  is  controlled  by  the 
user-specified  coupling-intensity  parameter  y  and  a  relaxation/ 
response  time,  t.  During  each  time  step  of  duration  At,  the 
particle-coordinate/cell  edge-length  scaling  factor,  p,  is  com¬ 
puted  as 

(1+TY(jP“P°))3  (EqB.l) 

where  P  is  the  instantaneous  pressure,  and  P0  is  the  target 
pressure.  The  instantaneous  pressure  is  obtained  from  the  viri- 
al  theorem  (Ref  44)  as 
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P  =  ^(K+W)  (EqB.2) 

where  K  and  V  are  respectively  the  instantaneous  kinetic 
energy  of  the  particles  residing  within  and  the  volume  of  the 
computational  cell,  and  W  is  the  instantaneous  virial  defined 
as  a  scalar  product  of  two  vectors,  one  defined  by  the  instan¬ 
taneous  coordinates  of  all  particles  in  the  cell,  and  the  other 
by  the  corresponding  force  components. 
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